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Abstract 

We have studied the properties of a classical A^-body system coupled to a bath containing Ng- 
body harmonic oscillators, employing an (Ns + Ng) model which is different from most of the 
existing models with N$ = 1. We have performed simulations for A^-oscillator systems, solving 
2(Ns + Ng) first-order differential equations with N$ — 1 — 10 and Ng ~ 10 — 1000, in order to 
calculate the time-dependent energy exchange between the system and the bath. The calculated 
energy in the system rapidly changes while its envelope has a much slower time dependence. 
Detailed calculations of the stationary energy distribution of the system fs(u) (u: an energy per 
particle in the system) have shown that its properties are mainly determined by N$ but weakly 
depend on Ng. The calculated fs{u) is analyzed with the use of the V and q-T distributions: the 
latter is derived with the superstatistical approach (SSA) and microcanonical approach (MCA) 
to the nonextensive statistics, where q stands for the entropic index. Based on analyses of our 
simulation results, a critical comparison is made between the SSA and MCA. Simulations have been 
performed also for the Ag-body ideal-gas system. The effect of the coupling between oscillators 
in the bath has been examined by additional (Ng + Ng) models which include baths consisting of 
coupled linear chains with periodic and fixed-end boundary conditions. 

PACS numbers: 05.40.-a, 05.70.-a, 05.10.Gg 
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I. INTRODUCTION 



The study on open systems is one of the important areas in classical and quantum statis- 
tics [I]. In the theory of open systems, the deterministic dynamics of particles in the system 
is replaced by the stochastic Langevin equation in the classical limit. The problem has been 
investigated with the use of various models in which a single particle (the system) is at- 
tached at the center (or edge) of a linear chain {2], 3|, or it is coupled to a bath consisting of 
a collection of harmonic oscillators (4]- 13]. Many studies have been made for open systems 
by using the Caldeira-Leggett (CL) model given by [4|-[6j 



CL 



EL 

2M 



V(Q) + J2 



n=l 



mtu, 



2m 2 



rt 



Qn nQ 



muj, 



(1) 



where M (m), P (p n ) and Q (q n ) denote the mass, momentum and coordinate of a particle 
in a system (bath), V(Q) stands for the potential in the system, u n the frequency of the 
nth oscillator in the iVg-body bath and c n the coupling constant between the system and 
bath. The CL model was originally introduced for infinite bath (Nb — > 00). In recent years, 
the CL model has been employed for a study of properties of a small system coupled to a 
finite bath joj-jl^]. A thermalization of a particle (the system) coupled to a finite bath has 
been investigated jjj, [lo| . It has been shown that a complete thermalization of the particle 
requires some conditions for relative ranges of oscillating frequencies in the system and bath 
|qI. Tlo|] - The specific heat of a single oscillator (the system) coupled to finite bath has been 



studied with the use of two different evaluation methods 
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121 ] . The energy exchange 
between particles in a rachet potential (the system) and finite bath (Nb = 1 —500) has been 



investigated 



13| 



Ford and Kac proposed the model given by 
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which is referred to as the FK model. The CL and FK models are formally equivalent {7] be- 
cause Eq. 02]) may be derived from Eq. ([T]) with c n = muj. However, the physical meanings 
of the coupling term in the CL and MK models are not the same. The CL model was initially 
introduced such that we take into account a linear coupling of — Q J2 n c nQn between system 
and bath Q], and then the counter term of c^Q 2 /mu^ was included for a compensation of 
the renormalization in the oscillating frequency by the introduced interaction. In contrast, 



the interaction term in Eq. (j2J) of the FK model clearly expresses the quadratic potential of 
springs between Q and q n . It is evident that the interaction term of the FK model in Eq. 
(T2J) preserves the translational invariance whereas that of the CL model in Eq. ([T]) does not 



in a strict sense 



14j-|l6| except for c n = mu^ for which the CL model reduces to the FK 



ationally invariant interaction in 
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16] . 



model as mentioned above. The importance of the trans 
the system plus bath models has been discussed in Refs. 

In existing models which have been proposed for open systems Q-[l3j], the number of 
particles in a systems is taken to be unity (N$ = 1) while a generic open system may 
contain any number of particles. It is necessary to develop an (N$ + Nb) model including 
finite iVg-body system (N$ > 1) coupled to iV^-body bath, with which we may investigate 
the properties of generic small systems. Extending the FK model, we will propose in this 
paper three types of (Ns + Nb) models (referred to as A, B and C). In the model A a bath 
consists of uncoupled oscillators, and in the models B and C baths contain coupled oscillators 
with the periodic and fixed-end boundary conditions, respectively. They are adopted for a 
study on effects of couplings in bath oscillators. 

In the last decade, many studies have been made for nonextensive statistics initially 
proposed by Tsallis [l3]-j2o|. In nonextensive systems, the probability distribution generally 
does not follow the Gaussian, but it is well described by the g-exponential distribution, 



p(u) cc e^=[l-(l- 9 )Ml /M , 



(3) 



where an inverse of the effective temperature /3q and the entropic index q are fitting param- 
eters, and the g-exponential function e x is defined by 



17|-|20| 



[1 + (1 - q)x 



(4) 



with [y} + = max(y, 0). In the limit of q — > 1.0, e x reduces to the exponential function e 
In a seminal paper [13], the g-exponential distribution was first derived by the maximum 



entropy method with the usejrf the so-called Tsallis entropy. Later superstatistical 
and microcanonical methods 
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24J have been proposed as alternative approaches to nonex- 



tensive statistics. Recent development has shown that small systems belong to nonextensive 
systems 20|]. Performing direct simulation (DS) for the proposed (Ns + Nb) model with 
the system containing independent Ns oscillators (the model A), we have calculated the 
stationary distribution of the system of fs(u) for the energy per particle u (= E s /N s , E s : 



the system energy). The calculated distribution is well described by the q-T distribution 
given by 



fs(u) oc u 



a—l—bu 



(5) 



where a, b and q are fitting parameters. It is easy to see that in the limit of q — > 1.0, the 
q-T distribution reduces to the conventional T distribution. As will be shown in Sec. Ill, 
superstatistical approach (SSA) [2l|, |22| and microcanonical approach (MCA) [23[-|29| lead 
to the equivalent expressions for fs(u) given by Eq. §5§ with a = N s and b = (3 N S , but 
with different expressions for the entropic index q: 



in the SSA, 
< 1.0 in the MCA. 



(6) 

<- i fi ;,, t iw> i\ if ' :\ 

{N B -1) 

The entropic index in the SSA is expressed in terms of a system parameter (Ns), while that 
in the MCA is expressed in terms of a bath parameter (Nb)- This difference is serious from 
the physical viewpoint of small open systems. The purpose of the present paper is twofold: 
to develop the (Ns + Nb) model in which an open system contains finite Ns particles, and 



21 



to investigate the validity of the stationary distribution functions derived in the SSA 
and MCA 23]- 29j. This is the first study on open systems with finite Ns (> 1) as far as 



we are aware of. 

The paper is organized as follows. In Sec. II, we propose the model A mentioned above, 
for which we perform DS of 2(Ns + Nb) differential equations for the As-oscillator system 
in order to calculate the time-dependent energy exchange between the system and bath. We 
present detailed calculations of fs(u), changing model parameters such as N s , N B , frequency 
distribution, mass of oscillators in the bath, and coupling strength between the system and 
bath. In Sec. Ill, we analyze the calculated fs(u) by using the T distribution [Eq. (136]) ] 
and the q-T distribution [Eq. ([5]) or ( |4"5]) ]. The former is derived based on the Boltzmann- 
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22j and MCA |23j-|29| of the 



Gibbs statistics and the latter is obtained with the SSA 
nonextensive statistics. DS has been made also for the system consisting of iVs-body ideal 
gases, whose results are compared to those of oscillators. We introduce the models B and C, 
whose DS for the oscillator systems will be reported. A comparison is made among Langevin 
equations derived in various models for open systems. The final Sec. IV is devoted to our 
conclusion. 
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II. ADOPTED (N s + N B ) MODEL 



A. A system with bath containing uncoupled oscillators 



We consider a system (H s ) and a bath (H B ) consisting of independent N s and N B one- 
dimensional oscillators, respectively, which are coupled by the interaction (Hi). We assume 
that the total Hamiltonian is given by 

H = H s + H B + Hi, 



with 
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(model A), 



(7) 

(8) 
(9) 
(10) 
(11) 



which is referred to as the model A. Here M (m) denotes the mass, Pk (p n ) the momentum, 
Qk (in) position of the oscillator, V(Qk) (v(q n )) the potential in the system (bath), c n k 
coupling constant, b n and u n spring constant and frequency in the bath, respectively, and 
f(t) an external force. A simple generalization of the FK model [Eq. (J2])] yields the model 
Hamiltonian given by Eq. ([7]) with Hs given by Eq. ([8]), H B = J2n=iPn/^ m an( ^ Hi = 
S/t=L S^=i( ma;2 /2)(Q , n — Qk) 2 - hi our model Hamiltonian, we have added v(q n ) in H B 
such that the Hamiltonian is symmetric with respect to an exchange of system -H- bath (for 
f(t) = 0) and such that we may discuss the coupled oscillators in baths (model B and C). 
Furthermore, we have included coupling Ck n in place of mu 2 in Hi of the generalized FK 
model in order to study the effect of system-bath couplings. We note that Hi in Eq. f fTUj) 
may be rewritten as 

N s /N B \ 1 N B / N s \ N s N B 

Hi = - ( Yl Ckn ) + 2 ( Ckn J 9 « ~ C knQkQn- (12) 

k=l \n=l J 71=1 \fc=l / k=l n=l 



Absorbing the first and second terms in Eq. (fl2j) to Hs and H B , respectively, we may 
regard the last term as the interaction. Such a model Hamiltonian with a linear coupling of 
— Sfc ^2 n c knQkQn corresponds to the generalized CL model for finite N s . 



From Eqs. (IT])- (fTT]) . we obtain 2(N$ + Nb) first-order differential equations, 
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(18) 



k=l 



prime (') and dot (•) denoting derivatives with respect to the argument and time, respectively. 
A formal solution of Eq. (118p for q n (t) is given by 

q n (t) = g n (0) cos u n t + ^" ^' sin u n t + / sin & n (t — t')Qe(1?) dt', 



with 



u n = — h> — = u n + y — . 

fe=i fc=i 



Substituting Eq. (jl9|) to Eq. (ITT)) , we obtain the Langevin equation given by 



Mft(i) = -V'(Q*)-M -J] / 
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(22) 
(23) 
(24) 



where denotes the additional interaction between k and £th particles in the system 
induced by couplings {ck n }, Jkt(t) the memory kernel and (k the stochastic force. 
If the equipartition relation is realized in initial values of q n (0) and q(0), 

(mu 2 n q n (0) 2 ) B = (mq n (0) 2 ) B = k B T, (25) 

we obtain the fluctuation-dissipation relation: 

(Ck(t)( k (t')) B = k B T lhk (t-t r ), (26) 

where (-) B stands for the average over variables in the bath. 

In the case of N B — > oo, summations in Eqs. ( 122|) - (|24|) are replaced by integrals. When 
the density of states (D(cu) = Ng 1 J2 n ^( u ~ w n)) is given by the Debye form: D{oj) oc to 2 
for < u < w B , the kernel becomes 

7 (t)oc f a5(t), (27) 

Tit 

which leads to the Markovian Langevin equation. 

In the case of N s = 1, we obtain £ and 7 in Eqs. f l2"2"|) and (|2"5|) where the subscripts k 
and £ are dropped (e.g., c kn = c n ), 

Nb / \ 

Mm = Em 1 -^? > ( 28 ) 



c; 



n=l 



1® = E K cos a; J. (29) 



The additional interaction vanishes (£ = 0) if we choose c n = mu 2 in Eq. (|28|) . 

In the case of iVg 7^ 1, however, it is impossible to choose {ck n } such that = is 
realized for all pairs of (k, t) in Eq. (122]) . Then is inevitably coupled to Qi for i 7^ k with 
the superexchange-type interaction of antiferromagnets: — J2 n c knC£ n /mtin m Eq. ( 1221) . 



B. Model calculations for oscillator systems 

It is easier to solve 2(N$ + N B ) first-order differential equations given by Eqs. ( fT3l) - 
( TT6l) than to solve the Ns Langevin equations given by Eqs. (|2T1) -(}!M1) although the latter 
provides us with clearer physical insight than the former. We have performed DS, solving 
the differential equations for the oscillator system with V{Qk) = MQlQl/2 in Eq. flH]) for 
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f(t) = 0, M = m = 1.0 and Qk = w„ = 1.0 otherwise noticed with the use of the fourth- 
order Runge-Kutta method with the time step of 0.01. In order to study the N$ and Nb 
dependences of various physical quantities, we have assumed the coupling given by 

Ckn = ivlk' (30) 

because the interaction term includes summations of J2k=i anc ^ ^2n=i * n ^q- i flOl - We have 
chosen Co = 10.0 (see Sec. II B 3. Effect of cq). It is noted that with our choice of Cfc n , the 
interaction contribution is finite even in the thermodynamical limit of Nb — > oo because 
the summation over n runs from 1 to Nb in Eq. ffTO]) . Although we have tried to adopt an 



alternative choice of Ck n given by 



3Q| 



(31) 



VNJNb~' 

qualitatively similar results have been obtained, as will be shown in Sec. Ill A. Initial 
conditions for Qfc(0), Qfc(0), q n (0) and q n (0) are given by random Gaussian variables with 
zero means and unit variances. Simulations have been performed for t = to 10000, results 
for t < 2000 being discarded for evaluations of stationary distributions. Results to be 
reported are averages over 10000 runs. 

We have assumed that the energies per particle u n {t) in the system (rj=S) and the bath 
(?7=B) are given by 



u s 



u B 



N B 

— Y 



2M 2 



2m 2 



(32) 



(33) 



neglecting a contribution from the interaction term Hj, which is valid for the weak interaction 



although a treatment of the finite interaction is ambiguous and controversial 
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12] . Figures 

Ufa) and (b) show the time dependence of u n for N$ = 1 and 10, respectively, with Nb = 1000 
of a single DS run. We note that although u v (t) rapidly oscillates, its envelope has much 
slower time dependence. Periods for rapid oscillations are about 0.95 and 2.22 for N$ = 1 and 
10, respectively: the latter value is larger than the former because of a larger renormalization 
effect due to couplings [the term in Eq. (1221)]. Magnitudes of time variations in us{t) are 
larger than those in «#(£) because N$ Nb- The width of variation in us{t) for N$ = 1 in 
Fig. QJa) is larger than that for N$ = 10 in Fig. QJb). Even when the energy of the system 
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is once decreased flowing into the bath, later it returns back to the system within the finite 

n 

time [31(. Then the dissipative energy transfer from the system to the bath or vice versa 
does not occur in a long time scale in Fig. [TJ This is in contrast with the result of Ref. 



131 ] which has reported a transition from non-dissipative to dissipative energy transfer at 
N B ~ 300 - 400 with N s = 1. 

In the following, we will show calculations of the stationary distributions of the system 
and bath, changing Ns, Nb, interaction strength (co), the distribution of u n and the ratio of 
m/M. Hereafter the argument u in the stationary distributions of fs(u) and /s(tt) expresses 
u = us and u = Ub, respectively. 

1. Effect of N s 

First we study the effect of Ns- Dashed, dotted, chain and solid curves in Fig. |2ta) show 
the stationary distribution of the system fs(u) for Ns = 1, 2, 5 and 10, respectively, with 
Nb = 100. fs(u) for Ns = 1 shows an exponential-like behavior while fs(u) for Ns > 1 
has a structure with a peak near the center of the stationary distribution of the bath 
Distributions of fs(u) for N s = 1, 2, 5 and 10 with N B = 100 are plotted by dashed, dotted, 
chain and solid curves, respectively, in Fig. 121b), which is nearly independent of Ns- More 
detailed discussion on the Ns dependence will be given in Sec. Ill A. 

2. Effect of N B 

Calculated distributions of fs(u) for Ns = 1 with Nb = 10, 1000 and 1000 are plotted by 
solid, dashed and chain curves, respectively, in Fig. Ela). Similar results of are shown 

in Fig. HHb). Profiles of fs(u) showing an exponential-like behavior are almost independent 
of Nb while those of change: its width becomes narrower for larger Nb- Solid, dashed 

and chain curves in Fig. g^a) [Fig. H^b)] show f s (u) [f B {u)} for N B = 10, 100 and 1000, 
respectively, with N$ = 10. Again fs{u) of Ns = 10 is nearly independent of Nb- In 
particular for Ns = Nb = 10, we obtain fs(u) = /b("u) because the system and bath are 
equivalent. for N$ = 10 in Fig. 4(b) is indistinguishable to that for N$ = 1 in Fig. 

3(b). 

3. Effect of c 

We change the coupling strength of c in Ck n = c /NsN B - Figure EI a), (b) and (c) show 
distributions of fs(u) and fs(u) for cq = 1.0, 10.0 and 100.0, respectively, with Ns = 1 and 
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N s = 10 for N B = 100. Results for c = 1.0 [Fig. GHa)] and c = 10.0 [Fig. EJb)] are almost 
identical. When Co is increased to 100.0, distribution of fs(u) becomes much wider than 
those in Figs. [5(a) and (b). At the same time, is modified by the stronger coupling. 

We have decided to adopt Co = 10.0 in our DS, related discussion being given in Sec. Ill A 
1. 

4- Effect of distributions of u n 

Although we have so far assumed u n = 1.0 in the bath, we will examine additional two 
types of distribution ranges for {u n }: uniform distributions in [0.5, 1.5] and [2.0,3.0] with a 
fixed flk = 1-0 in the system. Calculated fs(u) and f B {u) for u n G [0.5, 1.5] in Figs, E(b) 
are almost the same as those for u n = 1.0 in Figs. EJ^a). In Fig. [6(c) , where distribution of 
u n G [2.0,3.0] in the bath does not have an overlap with those of fl^ = 1.0 in the system, 
fs{u) is nearly the same as those in Figs. [6(a) and (b) in which the frequency ranges of the 
bath overlap those of the system. In contrast, in Fig. [6(c) is quite different from those 

in Figs. [6(a) and (b) as expected. Our results shown in Figs. [6(a), (b) and (c) suggest that 
fs(u) is not so sensitive to the position of frequency ranges of the bath relative to that of 
the system. This is in contrast with the result for Ns = 1 in Ref. jjj, which shows that for 
a thermalization of the system, the relative position between the oscillating frequency range 
of the system and that of the bath is very important. 

5. Effect ofm/M 

Finally we will change a value of m which has been so far assumed to be m = M = 1.0. 
Figures [7(a) and (b) show fs(u) for Ns = 1 and Ns = 10, respectively, with Nb = 100 for 
m/M = 1.0 (solid curves), 0.1 (dashed curves) and 0.01 (chain curves). With decreasing the 
ratio of m/M, the magnitude at small u of fs(u) is increased, by which that of fs(u) for 
Ns = 1 is slightly increased in Fig. [7J(a). However, the shape of fs(u) for Ns = 10 in Fig. 
Ulb) is almost unchanged with changing m/M. 

Before closing Sec. II, we point out that roles of the system and the bath are interchange- 
able in Figs. 2-7 because the Hamiltonian for the model A given by Eqs. ([7])- (llip has the 
system-B-bath symmetry for f(t) = 0. For example, /s(w) for Ns = 10 and Nb = 1 may be 
given by f B {u) in Fig. 3(b) for Ns = 1 and A?£ = 10. Figures 2-7 show that the properties 
of fs(u) are mainly determined by Ns, which is the main result of our study. 
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III. DISCUSSION 



A. Analysis of DS results for oscillator systems 

1. Boltzmann-Gibbs statistics 

We may theoretically evaluate the distribution of fs(u) as follows. First we calculate the 
distribution for a set of variables of {Qk, Vk) (Yk = Qk) with the Boltzmann-Gibbs statistics 
for the infinite bath characterized by the inverse temperature /3 (see Appendix A), 



f{Q,V) dQdV oc exp 



a ^(MV 2 MQ 2 Q 2 \ 

k=i v 7 



U^dQ k dV k , (34) 



oc Ep-h-^dEs, (35) 

where Eg denotes the energy in the system: E s = T,k=il MV k/ 2 + Mn 2 Q 2 /2}. From Eq. 
(13"5"j) . the distribution of the system for u (= E s /N s ) becomes 

fs(u) = ^u a - 1 e^ = g 1 {u), (36) 

with 

a = N S} b = N s (3 } (37) 

r(a) . , 

Z = ( 38 ) 

where g\{u) denotes the Y (or % 2 ) distribution, its subscript 1 being attached for later 
purpose. Mean (/i) and variance (a 2 ) of the T distribution are given by 

a 1 , s 

"-»-?• (39) 

a 1 , , . 

ff = F = j^- (40) 

Equation ( 139|) expresses the equipartition relation. The distribution fsiu) of the bath for u 
(= Eb/Nb) may be obtainable in a similar way where signifies the bath energy. 

Our DS in the preceding section has shown that the most influential parameter on the 
properties of the system is N$- We now pay our attention to the Ns dependence of calculated 
means (fj, v ) and root- mean-square (RMS) (a v ) of the system (r/ = S) and bath (r/ = B). Figure 
E] shows fi v and (7 CIS Sb function of Ns with Nb = 100 obtained by DS: filled (open) circles 



11 



denote [is {&s) of the system: filled (open) squares stand for fi B (a B ) of the bath. We obtain 
a s ) = (2.84, 3.77), (1.97, 2.12), (1.302, 0.755) and (1.097, 0.393) for N s = 1, 2, 5 and 10, 
respectively. With decreasing Ns from Ns = 10, fis and as are increased. In contrast, fi B 
and a B are almost independent of N$. An increase in /is with decreasing Ns is attributed 
to an increase in the effective frequency of the system given by [Eqs. (j22p and ( I30p ] 

1 



^fefc _ &kk 1" 



Co 



N s N 2 (N B u 2 + c /m)_ 

An increase in as with decreasing Ns is due to an increase in given by Eq. (|2"4"|) which 
is proportional to cq/Ns- When we adopt a smaller value of cq, these increases are reduced. 
For example, /is ((75) calculated with a smaller Co = 1.0 are plotted by filled (open) triangles 
in Fig. 8, which shows {fi s ,(Ts) = (1.12, 1.24), (1.02, 0.78), (0.97, 0.45) and (0.97, 0.31) for 
N s = 1, 2, 5 and 10, respectively. Related distributions for N s = 1 and 10 are plotted in 
Fig. 5(a). For this choice of Co = 1.0, the energy exchange between system and bath is 
considerably decreased. 

We have performed DS by using also an alternative choice of couplings of c kn = 
Cq/ V N S N B given by Eq. (l3Tj) . Filled and open diamonds in Fig. 8 show fi s and as, 
respectively, calculated with c' Q — 1/ vTO which is chosen such that Eq. fl3~T]) yields the same 
value of c kn = 0.01 as Eq. (ED} for N s = 10 and N B = 100. For N s = 1, 2, 5 and 10, 
we obtain (//<?, cr s ) = (1.61,2.10), (1.35,1.25), (1.18,0.58) and (1.097,0.393), respectively. 
With decreasing Ns, both fi s an d as are increased, which are qualitatively similar to those 
obtained with couplings given by Eq. ( 130|) . 

Next we examine the profiles of iVs-dependent fs(u)- By using the relation between 
parameters a and b in the T distribution with its average and variance given by Eqs. ( 1391) 
and (HDjl . we may determine a and b by a = ([i 2 /a 2 ) and b = fi/a 2 . With the use of 
the calculated /i s and a s , we obtain (a,b) = (0.565,0.199), (0.861,0.437), (2.98,2.28) and 
(7.78,7.11) for N s = 1, 2, 5 and 10, respectively. Unfortunately, these values of a are not 
in agreement with the theoretical value of a = Ns given by Eq. fl37|) . We have employed 
the T distribution given by Eq. ( 136]) with the parameters a and b determined above for our 
analysis of fs(u) shown in Fig. |2]^a). Dashed curves in Figs. [9]^a)-(d) show the calculated T 
distribution, while solid curves express DS results. We note in Figs. M^c) and (d) that the 
T distributions for Ns = 5 and 10 are ostensibly in good agreement with calculated fs{u) 
although the calculated a disagrees with the theoretical value of a (= iV^) as mentioned 
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above. Furthermore, an agreement becomes poor for results of smaller N$ = 1 and 2, 
whose analyses will be discussed with the use of the nonextensive statistics in the following 
sub-subsection. 



2. Nonextensive statistics 



a. Superstatistical approach 

A disagreement between theoretical results and DS ones might arise from a use of the 



Boltzmann-Gibbs statistics. We will analyze the calcu 



statistics 



17|-|20|. Wilk and Wlodarczk |2l| and Beck 



ated results by using the nonextensive 



22J have pointed out that the observed 



non-Gaussian distribution may be accounted for if we assume that the Gaussian distribution 
e _/3n is averaged over the T distribution of g((3) for fluctuating inverse temperature /3, 



p(u) oc e" 



3u g(f3)df3, 



(41) 



with 



( n 



n/2 



^n/2-l e -n/3/2/3 _ 



(42) 



r(n/2) V2/3 , 

Here n denotes the number of independent Gaussian Xi contributions to the \ 2 distribution 
of (3 = Y^i=i X? [22|, and /3q stands for the mean of (3: (3$ = ((3)p and variance is given by 



((3 2 )/3 — (3q = (2/n)/?Q. Equations f j4TT) and ( )42l) express the superstatistics whose concept 
may be understood such that complex nonextensive systems are in the nonequilibrium states 
with temporarily and spatially fluctuating inverse temperature. 

In order to more accurately account for our calculated fs(u), we employ the concept of 
the superstatistics. We assume that the T distribution f{u) given by Eq. (|36l) is averaged 
over the distribution g((3) given by Eq. (02]) with n = 2N$, 



N s -l-l3N s u 



fs(u) oc / (/-■■- ( 
'o 

U N S -1 

OC 



9(13) dfl, 



;i+/?o«) JVs ' 

With the normalization factor, fs(u) is expressed by the q-T distribution g q (u), 



(43) 
(44) 



fs(u) = - u a - l e~ bu 

Z <7 



9qW 



(45) 
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with 



1 



q 

a 
b 



1 - 
N s , 
N s (3 , 



(46) 
(47) 
(48) 



r(a)r(-i T - a ) 



[(q-i)b] a r(^- r ; 



r(a) 



i r(a ) r( 7 ^) 



for q > 1.0, 
for q = 1.0, 
for g < 1.0. 



(49) 



It is easy to see that in the limit of q — > 1.0, the q-T distribution g q {u) given by Eq. ( )45l) 

reduces to the V distribution gi(u) given by Eq. (1361) . Average and variance of the g-I" 
distribution are given by 

a 



a 



6[l-(g-l)(o + l)]' 

2 a ( 2 - <i) 



for q > 1.0, 



(50) 
(51) 



q b 2 [l - (q - l)(a + 2)][1 - (q - l)(a + l)] 2 
which reduce to /ii = a/b and cr 2 = a/b 2 for g = 1.0 in agreement with Eqs. (|39|) and (j4^j) . 
The g~r distribution g q (u) has a maximum at 

fa - 1) 



U = U r 



for a > 1.0. 



6[l-(g-l)(a-l)] 

The u dependence of g q {u) for typical parameters is shown in Appendix A (Fig. 15). 
b. Microcanonical approach I 



(52) 



Next we mention the MCA to the nonextensive statistics 23[]-29|]. We consider micro- 
canonical ensembles of particles with the energy E, which is divided into two subsystems 
1 and 2. A probability for subsystem 1 containing Aq particles to have energy E\ is given 
by Q,Q 

n 1 (E 1 )n 2 (E 2 ) 



/iVi (-^1 



(53) 



n 1+2 (E) ' 

where the structure function Q K (E) (k = 1, 2, 1 + 2) expresses the number of states with the 
energy E. We assume that fl K (E) is given by |23|, |27| |. 

Q K (E) = Km K E m -\ (54) 
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where K is a constant and m R the degrees of freedom of variables in subsystem k. Equation 
f )54|) is valid for ideal gases and harmonic oscillators with m K ^> 1. 

Interpreting subsystems 1 and 2 as a system and a bath, respectively, we apply the 
MCA to the oscillator system under consideration for which mj = N$ and mg = N B . For 
1 < N s < N B and E s < Eqs. QEJJ and (EU) yield 



/s(#s) oc £ 



1 - 



5 



1 - (1 - 



1/(1-4) 



Ns-l-PEs 
S e <5 ' 



(55) 

(56) 
(57) 



with 



q 



(N B -iy 



Nb 
E ' 



(58) 
(59) 



where we attach hats for quantities in the MCA to distinguish them from counterparts in 
the SSA. Equation ( 1571) is equivalent to the q-T distribution given by Eq. f T4"5l) if we read 
Es = Nsu and (3 = (3q. Similarly, we obtain the distribution defined by 23_|, \27 \ 

n B (E-E s ) 



Ps(E s ) 



oc e~ 



(60) 
(61) 



In the limit of N B — > 00, Eqs. (joYj) and (JbT|) reduce to 

f s (Es) oc - 1 e"^, 
Ps(E s ) oc e"^, 



(62) 
(63) 



with 



1 



(64) 



where the equipartition relation is employed for E B (3> Es). From Eqs. fl53l and fl60|) . a 
relation between f s (E s ) and p s (E s ) is given by 



/ S (E S ) = fi 5 (£ 5 )p 5 (£ 5 ). 
15 



(65) 



With increasing Eg, ps(Es) is decreased whereas Qs{Es) oc Eg S \ and then f s (E s ) has a 
maximum at E s = (N s - 1)//3[1 - (q - 1)(N S - 1)] for N s > 1. 



It should be noted that the g-exponential function adopted in Refs. |24j-[29] is defined 

by 

e x , = [1 + (q' - l)x) 1/iq '- l) for q' > 1, (66) 



which is different from that given by Eq. (J4j) proposed in Ref. [12] . The relation between q' 
and g is q' — 1 = 1 — q, with which Eq. ( 1581) becomes = 1 + 1/ (A^ — 1) (> 1.0). 

We have tried to apply the q-T distribution given by Eqs. (H5l) - (|49l) to an analysis of 
profiles of fs(u) in Fig. [9], but we could not obtain satisfactory results. Rather we have 
phenomenologically adopted the q-T distribution, choosing its parameters a, b and q such as 
to provide results in fairly good agreement with fs(u) in Fig. |9]with satisfying Eqs. ( |50l) and 
( 15"Tj) . Chain curves in Figs. [9^a) and (b) express g q {u) with (a,b,q) = (1.0,1.31,1.30) and 
(1.64, 1.36, 1.09), respectively, for Ns = 1 and 2, which have been tentatively determined by 
a cut and try method. It is note that fs(u) of DS is finite at u — 0.0 for N$ = 1, which 
requires a = 1.0. These chain curves are in better agreement with the calculated fs(u) than 
dashed curves expressing the T distribution. 

c. Microcanonical approach II 

We will derive the stationary distribution with the alternative MCA (MCA II). We again 
consider a collection of A^ particles with the energy E (= Meo) where eo denotes an appro- 
priate energy unit. A probability for its subsystem 1 containing A^ particles to have energy 
Ei (= Mxe ) is given by 

fm{Ml) ~ Wn(M) ' (67) 

with 

We apply Eq. (IBT)) to a system plus bath without using the condition: 1 C iVi <C JV, which 
is employed in the MCA I. We assume that M and Mi are real as given by 

N = N s + N B , Ni = N S , M = Es + Eb i Mi = ^, (69) 



16 



where E$ (Eb) denotes an energy in the system (bath). Then the probability for u 
Es/Ns) in the system is given by 

w Ns {M x ) w Nb (M - M x ) 



wn s +n b [M ) 



with 



ksT = \ = W B = "*> (71) 

where \Lb is the mean energy in the bath and wn{M) is given by Eq. fl68l) with a replacement 
of n! — > T(n + 1), T(x) being the T function. 

We have calculated fs(u) with the MCA II by using Eqs. floT?]) - fl7T]) . whose results with 
e = 1.0 (dashed curves), 0.1 (dotted curves) and 0.01 (chain curves) are shown in Figs. 
[TOTa)-(d). With decreasing 6q, results of MCA II are expected to approach the classical 
limit. Although a general trend is accounted for by MCA II calculations, their agreement 
with DS results is not so good. 



B. Comparison with ideal-gas systems 

Our model A given by Eqs. (I7|)- (ITT|) may be applied to ideal gases (the system) coupled 
to finite bath, for which we set V(Qj.) = 0. We have performed DS with the same model 
parameters (except for = 0) as in the case of oscillator systems mentioned in Sec. II. 
Solid curves in Figs. [TIT a), (b), (c) and (d) show calculated fs(u) of iV^-body ideal gases 
for Ns = 1, 2, 5 and 10, respectively, with iV^ = 100. For a comparison, we show by 
dotted curves, the corresponding results for oscillators having been plotted in Fig. [21(a). 
The energy distributions of bath, /b(w), for ideal-gas systems are almost the same as those 
for oscillator systems shown in Fig. [2(b). Comparing solid curves to dotted curves, we note 
that the distribution of fs(u) for ideal gases has larger magnitude at small u than that for 
oscillators. This yields the smaller average energy in ideal gases than that in oscillators, 
which is related with the fact the former has a smaller degree of freedom than the latter, as 
expressed in the equipartition relation. 

Figure 12 shows the Ns dependence of fi v and cr v = S,B: filled (open) circles show fis 
(as) and filled (open) squares denote jiB (°"b)- Although /i 5 in Fig. 12 has similar N s 
dependence to that in Fig. 8 for oscillator systems, magnitudes of the former are smaller 
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than the latter. The ratio of /jls(IG) / ns{OSC) approaches 0.5 with increasing N$, although 
the ratio is increased for Ns — > 1. 

We have analyzed calculated fs{u) in Figs. ITTT a)-(d) by using the T distribution given 
by Eq. d3S]) with 

a = ^, b = N s p, (72) 

which lead to 

" = T = W (73) 

°' = T? = 2fkp- (74) 

Equation ( 173|) expresses the equipartition relation of ideal gases. Simulations shown in Fig. 
11 yield Os,cr s ) = (2.59, 3.62), (1.50, 1.76), (0.716, 0.504) and (0.469, 0.222) for N s = 1, 
2, 5 and 10, respectively, from which we obtain (a, b) = (0.514 0.198), (0.733, 0.487), (2.02, 
2.82) and (4.47 9.54). Dashed curves in Figs. ITlT a)-(d) express the T distribution calculated 
with the use of a and b thus obtained. They are in good agreement with the DS results 
for Ns = 5 and 10, but not for Ns = 1 and 2. Chain curves express the q-T distributions 
obtained with (a, b, q) = (0.61, 0.38, 1.18) and (1.0, 0.85, 1.09) for N s = 1 and 2, respectively, 
which have been determined by a cut-and-try method. Results of the q-T distribution are 
in better agreement with DS than those of the T distribution. This situation is the same as 
in the case of oscillator systems as discussed in Sec. Ill A. 

C. Bath containing coupled oscillators 

In most of existing models for open systems jsj-jl^], baths are assumed to be consisting 
of uncoupled oscillators. In order to study the effect of couplings of oscillators in a bath, we 
consider the models B and C in which baths consist of coupled oscillators with the periodic 
and fixed-end boundary conditions, respectively. 

1. ModelB 

In the model B, we assume that the Hamiltonian is given by Eqs. (!7|)- (|T0|) with v(q n ), 

v (Qn) = ^(in - q n +i) 2 (modelB), (75) 
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under the periodic boundary condition: 

QN B +n = Qn, PN B +n = Pni 



(76) 



where b denotes the spring constant between neighboring sites in the bath and N B is assumed 
even without a loss of generality. 

Equations of motion for Q k and q n are given by 

Nb/2-1 

MQ k = -V'(Q k )- <Qk-qn)+f(t), (77) 



n=-N B /2 



N s 



mq n = -b(2q n - q n _ x - q n+l ) - ^ c (<?« - Qk)- 



(78) 



k=l 



By using a transformation mentioned in Appendix B, we obtain the Langevin equation 
for Qk(t) given by Eq. f[2"Tj) with 

c 2 N B 



Mf w = cN B 6 ki - 



moo, 



~,2 ' 



Jki(t) 
(k(t) 



c 2 N B 



~,2 



cvN B 



COS CUot, 

go(0) cos o)ot + sin a)ot 



with 



-.2 



cN s 



in 



(79) 
(80) 
(81) 

(82) 



2. Model C 

In the model C, we assume that the Hamiltonian is given by Eqs. (T71)- (|T0l) with v(q r , 

v (ln) = ^(<?n ~ <?n+i) 2 (modelC), 
under the fixed-end boundary condition given by 

% = (1n b +i =0, po = Pn b +i = 0. 



Equations of motion for Q k and q n are given by 



iV fl +l 



MQ k = -V'(Q k ) - <Qk-qn) + f(t), 



n=0 



N S 



mq n 



-b(2q n - q n _ x - q n+l ) - ^ c (<?n - Qk)- 



(83) 
(84) 

(85) 
(86) 



fc=i 
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By using a transformation mentioned in Appendix C, we obtain the Langevin equation 
given by Eq. (j2"Tl) with 



N 



Ck(t) 



c(N B + 2)5 M - 



B 11 

(rat 



s=l 



must 



N B 

£ 



, cosa; s i, 



N B 

' S ^ j ca s 

s=l 



g s (0) cosaXjt + saiQ K t 



Us 



(87) 
(88) 
(89) 



where Cj s and a s are expressed by 



c(iVs + 2) 

w. i 

m 



(90) 



2(N B + 1) 



cos 



7is (N B + 2)ns 
Y~ 2(N B + 1) 



, 7TS (AT B + 2)7TS 

<(,s, T + ^Tiy 



x cosec 



7T.S 



(91) 



2(N B + 1) 

DS calculations for models B and C have been performed for oscillator systems with the 
same parameters as in Sec. II in addition to b = 1.0. Dashed and solid curves in Fig. 
[151(a) [Fig. [T3fb)] show fs(u) [f B (u)] of the model B for Ng = 1 and 10, respectively, with 
N B = 100. Dashed and solid curves in Fig. [T3Ja) [Fig. [14|b)] show fs(u) [f B (u)] of the 
model C for Ns = 1 and 10, respectively, with iV^ = 100. Profiles of fs{u) and f B (u) of 
the model B in Fig. [13] are similar to those of the model C in Fig. [TH Comparing Figs. 
[TBI and [14] with Fig. [2J we note that couplings in oscillators of the bath have essentially no 
effects on the behavior of fs(u) of the system, although they have some effects on f B (u) as 
expected. 



D. Comparisons among various models 



Table 1 summarizes comparisons among elements of £kn,J,kn and (k in Langevin equations 
derived from various models for open systems including CL [5j] and MK models {7] and models 
A, B and C which are proposed in Sees. II and III. Additional interactions ^ n induced by 
introduced couplings between the system and bath remain finite in the models A, B and C 
although they vanish in the CL and MK models for Ns = 1. We note that functional forms 
of C,ke and (k in all the models are similar. This is the reason why properties of fs(u) and 
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TABLE I: Terms of 7^ and Cfc in the Langevin equation, MQk(t) = —V(Qk)—MY^£^kiQi(t)— 
Jo lki{t — t')Qe(t') dt' — 7fc£(*)Q£(0) + Ck(t), calculated by various models: 1) CL model [Eq. 
CQ]: 2) MK model [Eq. ©]: 3) the model A [Eq. ([TTjr j: 4) the model B [Eq. {75])]: 5) the model 
C [Eq. (]83j) ]. The CL and MK models are for Ns = 1 for which subscripts k,£ are dropped. 



model 




7fc<? 


Ck 


CL 1 ) (N s = 1) 





v ( 4 > 


COS U> n i 


En C ™ 


q n (0) C osu n t+ (2=M) 


sin cj„t 


MK 2 ) (Ns = 1) 





En mW n 


COS OJ n t 


En mW n 


<jr n (0) cos^ n t + (2=M 


J sinw n t 


A 3 ) 


Pt A, - c knC£„ 

Z^n \yknUkt mCj 2 




^ COS (it 


En C kn 


g„(0) cosw n i + 




B 4 ) 




V mw o J 


cVN^ 


■«,(0)ooB«2dt+(^) 


sin uj n t 


C 5) 


c(A B + 2)fo - E s Sf 


\^ mui'j 


COS CJ s t 




'&(0)cOBW 8 t+(^) 


sinw s t 



/b(m) in Figs. [21 [13] and [14] are similar. We note, however, that the kernel £^ of the model 
B is oscillating and not dissipative even for Nb — > 00, which arises from the translational 
symmetry in the bath. 



IV. CONCLUDING REMARKS 



It is worthwhile to make a comparison between the SSA and MCA, which lead to equiv- 
alent q-T distributions given by Eqs. (|4"5l) and ( 1571) . We should, however, note that the 
entropic index of fs(u) obtained in the SSA [Eq. (1461) ] is different from that derived in the 
MCA [Eq. (158]) ] as shown by Eq. ([6]): q in the SSA is expressed in terms of Ns and greater 
than unity, while q in the MCA is expressed in terms of Nb and less than unity. Our DS 
has shown that fs{u) depends on Ns in Fig. [2] or [9] while it is almost independent of Nb in 
Figs. [3] and HI which suggests that the entropic index of fs{u) depends mainly on Ns but 



only weakly on Nb- Furthermore, our phenomenological analyses show that t 
entropic indexes are greater than unity. These facts seem to support the SSA 



le deduced 
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22 1 but 



throw doubt on the MCA and its applications [23]- [29], although more detailed study is 
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necessary to draw a definite conclusion. 

To summarize, we have studied the properties of classical small systems coupled to finite 
bath, by employing the (N s + N B ) models A, B and C, in which iV^-body system is coupled 
to iV^-body bath. Simulations for oscillator and ideal-gas systems have shown the following: 

(i) the energy of the system oscillates rapidly although its envelope has much slower time 
dependence, 

(ii) the dissipation of the system energy is not observed in our DS with N$ ~ 1 — 10 and 
N B ~ 10 - 1000, 

(hi) the stationary energy distribution of the system fs{u) for N$ > 1 has a peak at about the 
average energy of the bath, although fs{u) for N$ = 1 has an exponential-like distribution 
decreasing monotonously with increasing u, 

(iv) calculated fs(u), whose properties depend mainly on Ns but only weakly on N B , may 
be phenomenologically described by the T or q-T distribution [Eq. (H5j) ]. and 

(v) the coupling among oscillators in the bath yields little effect in classical systems. 



The item (i) is consistent with a previous study for Ns = 1 in Ref. 131 ] . The item (ii) 



suggests that 
N B (> 1000) 



or the energy dissipation of system, we might need to adopt a much larger 



3l| ■ The thermalized state reported in Refs. 0, [lo| corresponds to our state 
for Ns = 1 with the exponential-like distribution, in agreement with the item (iii). The item 
(iv) is favorable to the SSA but not to the MCA although either of them cannot quantitatively 
explain the DS results. The item (v) is consistent with the classical specific heat of harmonic 
oscillators for which both Einstein and Debye models yield the same results. Our model A 
given by Eqs. (171)-(fTTT) is expected to have a wide applicability to classical small systems: for 
example, for studies on a system with various potentials V{Q) like the bi-stable potentials 
and on a work performed by time-dependent external force f(t) in Eq. ([8]). These subjects 
are left as our future study. 
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Appendix: A. q-x 2 and q-T distributions 



1. The q-x 2 distribution 

We will show that if n independent variables of {x{\ follow the g-Gaussian distribution, 
a variable defined by Y = Y^i=i x f follows the q-x 2 distribution with rank n defined by 

P(Y) = \ e- Y Y nl2 '\ (Al) 

where Z stands for the normalization factor. 

In order to derive Eq. (lAip . we first define a new variable of X 2 = Y^l=i x h ^ OT which we 
obtain 



p(x)dx oc e q i * J I dx 



i=l 

oc e-^X^dX, 

oc efY^'^Y-^dY, 



e -Y Y n/2-i dY , (A2) 



leading to the g-deformed x 2 distribution given by Eq. ( 1A1I) . 



It is noted that the factorization is not satisfied for the g-exponential function [17], l32| . 



e -M Y[e~ x , (A3) 

i 

except for q = 1 or n = 1. Then we cannot employ the method of the characteristic function 
by which the x 2 -function is conventionally derived from n independent Gaussian. 



2. The q-T distribution 

When generalizing n/2 in Eq. ( 1A1I) to a real number a, we obtain the q-T distribution, 

9 q (u) = ±r trt- 6 ", (A4) 

where Z q is given by Eq. ( )49l) . Some numerical examples of g q {u) are shown in Fig. [15j 
The q-T distribution for q > 1.0 has a larger magnitude than the T distri 
at large u because of the flat-tail properties of the g-exponential function 
g~r distribution for g < 1.0 has a compact structure because of cut-off properties of the 
g-exponential function with no magnitudes for u > 1/(1 — q)b. 
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Dution (q = 1.0) 
17| . In contrast, 



Appendix: B. Langevin equation in the model B 



We will explain a derivation of the Langevin equation in _t 
(17j)- ([TUj) . (1751) and ([75]) . By using the transformation given by 

ATs/2-l 



re model B given by Eqs. 



33|] 



(In 



Vn 



^ E 



i(2irns/N B )~ 



s=-N B /2 
JVb/2-1 



e i(2nns/N B )p^ 



s=-N B /2 



we obtain the diagonalized 



N B /2-l 



H 



B 



with 



s=-N B /2 



7T.S' 



E (t^i y ^ ■ -f < ■ 



9 Sin2 l^l (- = -^/2,..,^/2-l). 



Substituting Eqs. (JEU) and (|B2) to Eq. flTOJ lead to 



^5 N B /2~1 Ns 

fe=l s=-N B /2 k=l 



Then equations of motion become 



(Bl) 
(B2) 



(B3) 



(B4) 



(B5) 



MQ fc = -V'(Q k ) - cN B Q k + c^N B q Q + F(t), 

N S 

mq s = -mu% + c>/Nb~ ^ Qk^sO, 



with 



fc-i 



??? 



(B6) 
(B7) 



(B8) 



Note that the third term of Eq. (1B6I) and the second term of Eq. (1B7I) include only the 
s = component. Substituting a formal solution of q s to Eq. (1B6I) . we obtain the Langevin 
equation given by Eqs. fl2D and (T79]) - fl82]) . 
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Appendix: C. Langevin equation in the model C 



A derivation of the Langevin equation in the model C given by Eqs. (!7|)- (|T0|) . (!83|) and 
will be explained. A transformation given by js[ 3^ 

— Fi N B 



(In 



Pn 



N B + 



s=l 
N B 



nns 

sin | — r I q s , 



N B + 1 



t5> 



N B + 1^ \N B + 1 

s=l 



rms 

sin | — r I p s 



yields the diagonalized H B , 



N 



s=l 



2-2 



^ = 2.1^ + ^- 



with 



2 i' Ab \ ■ 2 
a;. = I — I sin 

m 



ITS 



2(N B + 1) 



[s = l,2,--N B ) 



From a transformation given by Eqs. ( 1C1I) and ( 1C2I) . we obtain Hi given by 
Hi = - 2^ ®k + ~ 2^ & ~ C Z^ ® k 1^ as<ls 



k=i 



k=l s=l 



with 



■N B +1 



N B + 



n=0 



sm 



Then equations of motion for and q s become 



ixns 
N B + 1 



MQ k 



-V'{Q k ) - c{N B + 2)Q k + cJ2 a S q S + F{t), 



s=l 



N S 



mq s = -mu 2 q s + ca s 2J Qk, 



with 



(CI) 
(C2) 



(C3) 



(C4) 



(C5) 



(C6) 



(C7) 
(C8) 



mul = moj 2 s + c(N s + 2), 



(C9) 



Substituting a formal solution of q s to Eq. (1C7j) . we obtain the Langevin equation given by 
Eqs. (J2U) and (IH7D-(|89D. 
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The s dependence of a s given by Eq. (1C6j) or Eq. (I9T1) is plotted in Fig. [161 showing the 
zig-toothed structure whose magnitude decreases rapidly with increasing s. 
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FIG. 1: (Color online) Time dependences of us(t) and ue(f) for (a) N$ = 1 and (b) N$ = 10 with 
Nb = 1000 (a single DS run), inset showing enlarged plots of us(t) for t = to 60. 

FIG. 2: (Color online) (a) Stationary distributions of (a) fs(u) and (b) /b(u) with iVe = 100 for 
various Ng: N$ = 1 (dashed curves), 2 (dotted curves), 5 (chain curves) and 10 (solid curves). 

FIG. 4: (Color online) Stationary distributions of (a) fs(u) and (b) /b(u) with N$ = 10 for 
various Nb'- Nb = 10 (solid curves), 100 (dashed curves) and 1000 (chain curves). 

FIG. 5: (Color online) Stationary distributions of fs(u) and for (a) Co = 1.0, (b) 10.0 and 

(c) 100.0 in Ckn = c /NsNb: chain (solid) curve denotes fs for N$ = 1 (Ns = 10), and dotted 
(dashed) curve expresses fs for Ng = 1 (Ns = 10) with Nb = 100, Jb(u) being divided by a factor 
of two. 

FIG. 6: (Color online) Stationary distributions of fs(u) (solid curves) and /s(m) (dashed curves) 
for (a) oj n = 1.0, (b) u n G [0.5,1.5] and (c) uj n G [2.0,3.0] with N s = 10 and N B = 100, f B (u) 
being divided by a factor of two. 

FIG. 7: (Color online) Stationary distributions of fs(u) for (a) Ns = 1 and (b) Ns = 10 with 
Nb = 100 for various m/M: m/M = 1.0 (solid curves), 0.1 (dashed curves) and 0.01 (chain curves). 

FIG. 3: (Color online) Stationary distributions of (a) fs{u) and (b) /b{u) with Ns = 1 for various 
Nb- Nb = 10 (solid curves), 100 (dashed curves) and 1000 (chain curves). 
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FIG. 8: (Color online) N$ dependences of (i^ and of systems (rj= S) and baths (rj= B) with 
Nb = 100: filled (open) circles show fig (&s)> and filled (open) squares \ib (&b) with cq = 10.0: 
filled (open) triangles express fM$ (°\s) calculated with cq = 1.0: filled (open) diamonds denote us 
(as) calculated with the coupling given by Cfc n = 1.0/\/WNsNb [Eq. ([3"T]) ] (see text). 

FIG. 9: (Color online) The u dependence of f s (u) for (a) N s = 1, (b) Afe = 2, (c) iV s = 5 
and (d) Ns = 10 with Nb = 100 obtained by our direct simulation (DS: solid curves) and the T 
distribution (r) given by Eq. (|36[) (dashed curved). Chain curves in (a) and (b) express the q-T 
distribution (q-T) given by Eq. ([35]) (see text). 

FIG. 10: (Color online) The u dependence of f s (u) for (a) N s = 1, (b) N s = 2, (c) = 5 and 
(d) N s = 10 with Nb = 100 obtained by the MCA II with e = 1.0 (dashed curves), 0.1 (dotted 
curves) and 0.01 (chain curves) [Eqs. ([69]) - ([7T]) ]. solid curves expressing DS results. 

FIG. 11: (Color online) The u dependence of fs(u) of ideal-gas systems for (a) Ns = 1, (b) 
N s = 2, (c) N s = 5 and (d) N s = 10 with N B = 100: DS (IG: solid curves), the T distribution 
(r (IG): dashed curves) and the q-T distribution (q-T (IG): chain curves). For a comparison, DS 
results for oscillator system (OSC) are plotted by dotted curves (see text). 

FIG. 12: (Color online) Ns dependences of ^ and a v of ideal gas systems (ry= S: circles) and 
baths (77= B: squares) with iV^ = 100: filled and open marks denote mean and RMS, respectively. 

FIG. 13: (Color online) Stationary distributions of (a) fs(u) and (b) /b(u) of the model B for 
Ns = 1 (dashed curves) and 10 (solid curves) with A^^ = 100. 
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FIG. 14: (Color online) Stationary distributions of (a) fs(u) and (b) /b(u) of the model C for 
Ns = 1 (dashed curves) and 10 (solid curves) with Nb = 100. 



FIG. 15: (Color online) The q-T distribution g q (u) [Eq. (j%4j) ] for (a) a = 0.5, (b) 1.0, (c) 1.5, 
(d) 2.0, (e) 3.0 and (f) 4.0 with b = 1.0: q = 0.9 (chain curves), 1.0 (dashed curves) and 1.1 (solid 
curves) . 



FIG. 16: (Color online) The s dependence of a s for N = 10 (dashed curve), 20 (chain curve), 50 
(dotted curve) and 100 (solid curve) [Eq. ((C6]) or (|9"T]) ]. 
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